// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_REDUX_H
#define EIGEN_REDUX_H

namespace Eigen {

namespace internal {

// TODO
//  * implement other kind of vectorization
//  * factorize code

/***************************************************************************
 * Part 1 : the logic deciding a strategy for vectorization and unrolling
 ***************************************************************************/

template<typename Func, typename Evaluator>
struct redux_traits
{
  public:
	typedef typename find_best_packet<typename Evaluator::Scalar, Evaluator::SizeAtCompileTime>::type PacketType;
	enum
	{
		PacketSize = unpacket_traits<PacketType>::size,
		InnerMaxSize = int(Evaluator::IsRowMajor) ? Evaluator::MaxColsAtCompileTime : Evaluator::MaxRowsAtCompileTime,
		OuterMaxSize = int(Evaluator::IsRowMajor) ? Evaluator::MaxRowsAtCompileTime : Evaluator::MaxColsAtCompileTime,
		SliceVectorizedWork = int(InnerMaxSize) == Dynamic	 ? Dynamic
							  : int(OuterMaxSize) == Dynamic ? (int(InnerMaxSize) >= int(PacketSize) ? Dynamic : 0)
															 : (int(InnerMaxSize) / int(PacketSize)) * int(OuterMaxSize)
	};

	enum
	{
		MightVectorize = (int(Evaluator::Flags) & ActualPacketAccessBit) && (functor_traits<Func>::PacketAccess),
		MayLinearVectorize = bool(MightVectorize) && (int(Evaluator::Flags) & LinearAccessBit),
		MaySliceVectorize =
			bool(MightVectorize) && (int(SliceVectorizedWork) == Dynamic || int(SliceVectorizedWork) >= 3)
	};

  public:
	enum
	{
		Traversal = int(MayLinearVectorize)	 ? int(LinearVectorizedTraversal)
					: int(MaySliceVectorize) ? int(SliceVectorizedTraversal)
											 : int(DefaultTraversal)
	};

  public:
	enum
	{
		Cost = Evaluator::SizeAtCompileTime == Dynamic
				   ? HugeCost
				   : int(Evaluator::SizeAtCompileTime) * int(Evaluator::CoeffReadCost) +
						 (Evaluator::SizeAtCompileTime - 1) * functor_traits<Func>::Cost,
		UnrollingLimit = EIGEN_UNROLLING_LIMIT * (int(Traversal) == int(DefaultTraversal) ? 1 : int(PacketSize))
	};

  public:
	enum
	{
		Unrolling = Cost <= UnrollingLimit ? CompleteUnrolling : NoUnrolling
	};

#ifdef EIGEN_DEBUG_ASSIGN
	static void debug()
	{
		std::cerr << "Xpr: " << typeid(typename Evaluator::XprType).name() << std::endl;
		std::cerr.setf(std::ios::hex, std::ios::basefield);
		EIGEN_DEBUG_VAR(Evaluator::Flags)
		std::cerr.unsetf(std::ios::hex);
		EIGEN_DEBUG_VAR(InnerMaxSize)
		EIGEN_DEBUG_VAR(OuterMaxSize)
		EIGEN_DEBUG_VAR(SliceVectorizedWork)
		EIGEN_DEBUG_VAR(PacketSize)
		EIGEN_DEBUG_VAR(MightVectorize)
		EIGEN_DEBUG_VAR(MayLinearVectorize)
		EIGEN_DEBUG_VAR(MaySliceVectorize)
		std::cerr << "Traversal"
				  << " = " << Traversal << " (" << demangle_traversal(Traversal) << ")" << std::endl;
		EIGEN_DEBUG_VAR(UnrollingLimit)
		std::cerr << "Unrolling"
				  << " = " << Unrolling << " (" << demangle_unrolling(Unrolling) << ")" << std::endl;
		std::cerr << std::endl;
	}
#endif
};

/***************************************************************************
 * Part 2 : unrollers
 ***************************************************************************/

/*** no vectorization ***/

template<typename Func, typename Evaluator, int Start, int Length>
struct redux_novec_unroller
{
	enum
	{
		HalfLength = Length / 2
	};

	typedef typename Evaluator::Scalar Scalar;

	EIGEN_DEVICE_FUNC
	static EIGEN_STRONG_INLINE Scalar run(const Evaluator& eval, const Func& func)
	{
		return func(redux_novec_unroller<Func, Evaluator, Start, HalfLength>::run(eval, func),
					redux_novec_unroller<Func, Evaluator, Start + HalfLength, Length - HalfLength>::run(eval, func));
	}
};

template<typename Func, typename Evaluator, int Start>
struct redux_novec_unroller<Func, Evaluator, Start, 1>
{
	enum
	{
		outer = Start / Evaluator::InnerSizeAtCompileTime,
		inner = Start % Evaluator::InnerSizeAtCompileTime
	};

	typedef typename Evaluator::Scalar Scalar;

	EIGEN_DEVICE_FUNC
	static EIGEN_STRONG_INLINE Scalar run(const Evaluator& eval, const Func&)
	{
		return eval.coeffByOuterInner(outer, inner);
	}
};

// This is actually dead code and will never be called. It is required
// to prevent false warnings regarding failed inlining though
// for 0 length run() will never be called at all.
template<typename Func, typename Evaluator, int Start>
struct redux_novec_unroller<Func, Evaluator, Start, 0>
{
	typedef typename Evaluator::Scalar Scalar;
	EIGEN_DEVICE_FUNC
	static EIGEN_STRONG_INLINE Scalar run(const Evaluator&, const Func&) { return Scalar(); }
};

/*** vectorization ***/

template<typename Func, typename Evaluator, int Start, int Length>
struct redux_vec_unroller
{
	template<typename PacketType>
	EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE PacketType run(const Evaluator& eval, const Func& func)
	{
		enum
		{
			PacketSize = unpacket_traits<PacketType>::size,
			HalfLength = Length / 2
		};

		return func.packetOp(
			redux_vec_unroller<Func, Evaluator, Start, HalfLength>::template run<PacketType>(eval, func),
			redux_vec_unroller<Func, Evaluator, Start + HalfLength, Length - HalfLength>::template run<PacketType>(
				eval, func));
	}
};

template<typename Func, typename Evaluator, int Start>
struct redux_vec_unroller<Func, Evaluator, Start, 1>
{
	template<typename PacketType>
	EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE PacketType run(const Evaluator& eval, const Func&)
	{
		enum
		{
			PacketSize = unpacket_traits<PacketType>::size,
			index = Start * PacketSize,
			outer = index / int(Evaluator::InnerSizeAtCompileTime),
			inner = index % int(Evaluator::InnerSizeAtCompileTime),
			alignment = Evaluator::Alignment
		};
		return eval.template packetByOuterInner<alignment, PacketType>(outer, inner);
	}
};

/***************************************************************************
 * Part 3 : implementation of all cases
 ***************************************************************************/

template<typename Func,
		 typename Evaluator,
		 int Traversal = redux_traits<Func, Evaluator>::Traversal,
		 int Unrolling = redux_traits<Func, Evaluator>::Unrolling>
struct redux_impl;

template<typename Func, typename Evaluator>
struct redux_impl<Func, Evaluator, DefaultTraversal, NoUnrolling>
{
	typedef typename Evaluator::Scalar Scalar;

	template<typename XprType>
	EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Scalar run(const Evaluator& eval, const Func& func, const XprType& xpr)
	{
		eigen_assert(xpr.rows() > 0 && xpr.cols() > 0 && "you are using an empty matrix");
		Scalar res;
		res = eval.coeffByOuterInner(0, 0);
		for (Index i = 1; i < xpr.innerSize(); ++i)
			res = func(res, eval.coeffByOuterInner(0, i));
		for (Index i = 1; i < xpr.outerSize(); ++i)
			for (Index j = 0; j < xpr.innerSize(); ++j)
				res = func(res, eval.coeffByOuterInner(i, j));
		return res;
	}
};

template<typename Func, typename Evaluator>
struct redux_impl<Func, Evaluator, DefaultTraversal, CompleteUnrolling>
	: redux_novec_unroller<Func, Evaluator, 0, Evaluator::SizeAtCompileTime>
{
	typedef redux_novec_unroller<Func, Evaluator, 0, Evaluator::SizeAtCompileTime> Base;
	typedef typename Evaluator::Scalar Scalar;
	template<typename XprType>
	EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Scalar run(const Evaluator& eval,
															const Func& func,
															const XprType& /*xpr*/)
	{
		return Base::run(eval, func);
	}
};

template<typename Func, typename Evaluator>
struct redux_impl<Func, Evaluator, LinearVectorizedTraversal, NoUnrolling>
{
	typedef typename Evaluator::Scalar Scalar;
	typedef typename redux_traits<Func, Evaluator>::PacketType PacketScalar;

	template<typename XprType>
	static Scalar run(const Evaluator& eval, const Func& func, const XprType& xpr)
	{
		const Index size = xpr.size();

		const Index packetSize = redux_traits<Func, Evaluator>::PacketSize;
		const int packetAlignment = unpacket_traits<PacketScalar>::alignment;
		enum
		{
			alignment0 = (bool(Evaluator::Flags & DirectAccessBit) && bool(packet_traits<Scalar>::AlignedOnScalar))
							 ? int(packetAlignment)
							 : int(Unaligned),
			alignment = EIGEN_PLAIN_ENUM_MAX(alignment0, Evaluator::Alignment)
		};
		const Index alignedStart = internal::first_default_aligned(xpr);
		const Index alignedSize2 = ((size - alignedStart) / (2 * packetSize)) * (2 * packetSize);
		const Index alignedSize = ((size - alignedStart) / (packetSize)) * (packetSize);
		const Index alignedEnd2 = alignedStart + alignedSize2;
		const Index alignedEnd = alignedStart + alignedSize;
		Scalar res;
		if (alignedSize) {
			PacketScalar packet_res0 = eval.template packet<alignment, PacketScalar>(alignedStart);
			if (alignedSize > packetSize) // we have at least two packets to partly unroll the loop
			{
				PacketScalar packet_res1 = eval.template packet<alignment, PacketScalar>(alignedStart + packetSize);
				for (Index index = alignedStart + 2 * packetSize; index < alignedEnd2; index += 2 * packetSize) {
					packet_res0 = func.packetOp(packet_res0, eval.template packet<alignment, PacketScalar>(index));
					packet_res1 =
						func.packetOp(packet_res1, eval.template packet<alignment, PacketScalar>(index + packetSize));
				}

				packet_res0 = func.packetOp(packet_res0, packet_res1);
				if (alignedEnd > alignedEnd2)
					packet_res0 =
						func.packetOp(packet_res0, eval.template packet<alignment, PacketScalar>(alignedEnd2));
			}
			res = func.predux(packet_res0);

			for (Index index = 0; index < alignedStart; ++index)
				res = func(res, eval.coeff(index));

			for (Index index = alignedEnd; index < size; ++index)
				res = func(res, eval.coeff(index));
		} else // too small to vectorize anything.
			   // since this is dynamic-size hence inefficient anyway for such small sizes, don't try to optimize.
		{
			res = eval.coeff(0);
			for (Index index = 1; index < size; ++index)
				res = func(res, eval.coeff(index));
		}

		return res;
	}
};

// NOTE: for SliceVectorizedTraversal we simply bypass unrolling
template<typename Func, typename Evaluator, int Unrolling>
struct redux_impl<Func, Evaluator, SliceVectorizedTraversal, Unrolling>
{
	typedef typename Evaluator::Scalar Scalar;
	typedef typename redux_traits<Func, Evaluator>::PacketType PacketType;

	template<typename XprType>
	EIGEN_DEVICE_FUNC static Scalar run(const Evaluator& eval, const Func& func, const XprType& xpr)
	{
		eigen_assert(xpr.rows() > 0 && xpr.cols() > 0 && "you are using an empty matrix");
		const Index innerSize = xpr.innerSize();
		const Index outerSize = xpr.outerSize();
		enum
		{
			packetSize = redux_traits<Func, Evaluator>::PacketSize
		};
		const Index packetedInnerSize = ((innerSize) / packetSize) * packetSize;
		Scalar res;
		if (packetedInnerSize) {
			PacketType packet_res = eval.template packet<Unaligned, PacketType>(0, 0);
			for (Index j = 0; j < outerSize; ++j)
				for (Index i = (j == 0 ? packetSize : 0); i < packetedInnerSize; i += Index(packetSize))
					packet_res =
						func.packetOp(packet_res, eval.template packetByOuterInner<Unaligned, PacketType>(j, i));

			res = func.predux(packet_res);
			for (Index j = 0; j < outerSize; ++j)
				for (Index i = packetedInnerSize; i < innerSize; ++i)
					res = func(res, eval.coeffByOuterInner(j, i));
		} else // too small to vectorize anything.
			   // since this is dynamic-size hence inefficient anyway for such small sizes, don't try to optimize.
		{
			res = redux_impl<Func, Evaluator, DefaultTraversal, NoUnrolling>::run(eval, func, xpr);
		}

		return res;
	}
};

template<typename Func, typename Evaluator>
struct redux_impl<Func, Evaluator, LinearVectorizedTraversal, CompleteUnrolling>
{
	typedef typename Evaluator::Scalar Scalar;

	typedef typename redux_traits<Func, Evaluator>::PacketType PacketType;
	enum
	{
		PacketSize = redux_traits<Func, Evaluator>::PacketSize,
		Size = Evaluator::SizeAtCompileTime,
		VectorizedSize = (int(Size) / int(PacketSize)) * int(PacketSize)
	};

	template<typename XprType>
	EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Scalar run(const Evaluator& eval, const Func& func, const XprType& xpr)
	{
		EIGEN_ONLY_USED_FOR_DEBUG(xpr)
		eigen_assert(xpr.rows() > 0 && xpr.cols() > 0 && "you are using an empty matrix");
		if (VectorizedSize > 0) {
			Scalar res = func.predux(
				redux_vec_unroller<Func, Evaluator, 0, Size / PacketSize>::template run<PacketType>(eval, func));
			if (VectorizedSize != Size)
				res = func(
					res, redux_novec_unroller<Func, Evaluator, VectorizedSize, Size - VectorizedSize>::run(eval, func));
			return res;
		} else {
			return redux_novec_unroller<Func, Evaluator, 0, Size>::run(eval, func);
		}
	}
};

// evaluator adaptor
template<typename _XprType>
class redux_evaluator : public internal::evaluator<_XprType>
{
	typedef internal::evaluator<_XprType> Base;

  public:
	typedef _XprType XprType;
	EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE explicit redux_evaluator(const XprType& xpr)
		: Base(xpr)
	{
	}

	typedef typename XprType::Scalar Scalar;
	typedef typename XprType::CoeffReturnType CoeffReturnType;
	typedef typename XprType::PacketScalar PacketScalar;

	enum
	{
		MaxRowsAtCompileTime = XprType::MaxRowsAtCompileTime,
		MaxColsAtCompileTime = XprType::MaxColsAtCompileTime,
		// TODO we should not remove DirectAccessBit and rather find an elegant way to query the alignment offset at
		// runtime from the evaluator
		Flags = Base::Flags & ~DirectAccessBit,
		IsRowMajor = XprType::IsRowMajor,
		SizeAtCompileTime = XprType::SizeAtCompileTime,
		InnerSizeAtCompileTime = XprType::InnerSizeAtCompileTime
	};

	EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffReturnType coeffByOuterInner(Index outer, Index inner) const
	{
		return Base::coeff(IsRowMajor ? outer : inner, IsRowMajor ? inner : outer);
	}

	template<int LoadMode, typename PacketType>
	EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE PacketType packetByOuterInner(Index outer, Index inner) const
	{
		return Base::template packet<LoadMode, PacketType>(IsRowMajor ? outer : inner, IsRowMajor ? inner : outer);
	}
};

} // end namespace internal

/***************************************************************************
 * Part 4 : public API
 ***************************************************************************/

/** \returns the result of a full redux operation on the whole matrix or vector using \a func
 *
 * The template parameter \a BinaryOp is the type of the functor \a func which must be
 * an associative operator. Both current C++98 and C++11 functor styles are handled.
 *
 * \warning the matrix must be not empty, otherwise an assertion is triggered.
 *
 * \sa DenseBase::sum(), DenseBase::minCoeff(), DenseBase::maxCoeff(), MatrixBase::colwise(), MatrixBase::rowwise()
 */
template<typename Derived>
template<typename Func>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
DenseBase<Derived>::redux(const Func& func) const
{
	eigen_assert(this->rows() > 0 && this->cols() > 0 && "you are using an empty matrix");

	typedef typename internal::redux_evaluator<Derived> ThisEvaluator;
	ThisEvaluator thisEval(derived());

	// The initial expression is passed to the reducer as an additional argument instead of
	// passing it as a member of redux_evaluator to help
	return internal::redux_impl<Func, ThisEvaluator>::run(thisEval, func, derived());
}

/** \returns the minimum of all coefficients of \c *this.
 * In case \c *this contains NaN, NaNPropagation determines the behavior:
 *   NaNPropagation == PropagateFast : undefined
 *   NaNPropagation == PropagateNaN : result is NaN
 *   NaNPropagation == PropagateNumbers : result is minimum of elements that are not NaN
 * \warning the matrix must be not empty, otherwise an assertion is triggered.
 */
template<typename Derived>
template<int NaNPropagation>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
DenseBase<Derived>::minCoeff() const
{
	return derived().redux(Eigen::internal::scalar_min_op<Scalar, Scalar, NaNPropagation>());
}

/** \returns the maximum of all coefficients of \c *this.
 * In case \c *this contains NaN, NaNPropagation determines the behavior:
 *   NaNPropagation == PropagateFast : undefined
 *   NaNPropagation == PropagateNaN : result is NaN
 *   NaNPropagation == PropagateNumbers : result is maximum of elements that are not NaN
 * \warning the matrix must be not empty, otherwise an assertion is triggered.
 */
template<typename Derived>
template<int NaNPropagation>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
DenseBase<Derived>::maxCoeff() const
{
	return derived().redux(Eigen::internal::scalar_max_op<Scalar, Scalar, NaNPropagation>());
}

/** \returns the sum of all coefficients of \c *this
 *
 * If \c *this is empty, then the value 0 is returned.
 *
 * \sa trace(), prod(), mean()
 */
template<typename Derived>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
DenseBase<Derived>::sum() const
{
	if (SizeAtCompileTime == 0 || (SizeAtCompileTime == Dynamic && size() == 0))
		return Scalar(0);
	return derived().redux(Eigen::internal::scalar_sum_op<Scalar, Scalar>());
}

/** \returns the mean of all coefficients of *this
 *
 * \sa trace(), prod(), sum()
 */
template<typename Derived>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
DenseBase<Derived>::mean() const
{
#ifdef __INTEL_COMPILER
#pragma warning push
#pragma warning(disable : 2259)
#endif
	return Scalar(derived().redux(Eigen::internal::scalar_sum_op<Scalar, Scalar>())) / Scalar(this->size());
#ifdef __INTEL_COMPILER
#pragma warning pop
#endif
}

/** \returns the product of all coefficients of *this
 *
 * Example: \include MatrixBase_prod.cpp
 * Output: \verbinclude MatrixBase_prod.out
 *
 * \sa sum(), mean(), trace()
 */
template<typename Derived>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
DenseBase<Derived>::prod() const
{
	if (SizeAtCompileTime == 0 || (SizeAtCompileTime == Dynamic && size() == 0))
		return Scalar(1);
	return derived().redux(Eigen::internal::scalar_product_op<Scalar>());
}

/** \returns the trace of \c *this, i.e. the sum of the coefficients on the main diagonal.
 *
 * \c *this can be any matrix, not necessarily square.
 *
 * \sa diagonal(), sum()
 */
template<typename Derived>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
MatrixBase<Derived>::trace() const
{
	return derived().diagonal().sum();
}

} // end namespace Eigen

#endif // EIGEN_REDUX_H
